
############################################
##### Simulate Effects from Main Model #####
############################################

####################
##### Figure 1 #####
####################

library(foreign)
library(lme4)
library(scales)
library(ggplot2)
library(memisc)
library(Zelig)
library(MASS)

# real data
setwd("~")
admin <- read.dta("~services_admin_tscs.dta")
admin <- admin[admin$region_x=="Sub-Saharan Africa",]

# coefs.dta and covariance.txt produced in Main_Results_Replication_R.R

# parameter estimates
parm <- read.dta("coefs.dta")

means <- parm$estimate

#variance-covariance matrix
var <- read.table("covariance.txt")


#data means of control variables
varmeans <- c(mean(admin$lpop_l,na.rm=T),mean(admin$wdi_urban_l,na.rm=T),mean(admin$lgdppc_l,na.rm=T),
              mean(admin$conflict_l,na.rm=T),mean(admin$dpi_state_l,na.rm=T),mean(admin$p_polity2_l,na.rm=T),
              mean(admin$loilpc_l,na.rm=T),mean(admin$aid_pc_l,na.rm=T),rep(0,33),1,1)

#draws from a multivariate normal
draws <- mvrnorm(n=1000, means, var)
max <- as.numeric(quantile(admin$adminpc,na.rm=T,probs=c(1)))
max <- 60

adminvec <- seq(min(admin$adminpc,na.rm=T),max, by=((max)-min(admin$adminpc,na.rm=T))/(500-1))

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(adminvec)){
  out <- draws[,43] + varmeans[1]*draws[,1]+ varmeans[2]*draws[,2]+varmeans[3]*draws[,3]+varmeans[4]*draws[,4]+varmeans[5]*draws[,5]+varmeans[6]*draws[,6]+varmeans[7]*draws[,7]+varmeans[8]*draws[,8]+adminvec[i]*draws[,9]+ (adminvec[i]^2)*draws[,10]+varmeans[41]*draws[,41]
  ev <- out
  meanhelp <- means[43] + varmeans[1]*means[1]+ varmeans[2]*means[2]+varmeans[3]*means[3]+varmeans[4]*means[4]+varmeans[5]*means[5]+varmeans[6]*means[6]+varmeans[7]*means[7]+varmeans[8]*means[8]+adminvec[i]*means[9]+ (adminvec[i]^2)*means[10] + varmeans[41]*means[41]
  mean <- meanhelp
  lower <- quantile(ev,0.05)
  higher <- quantile(ev,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

plots <- data.frame(adminvec,meanvec)
plotadmin <- data.frame(admin$adminpc[admin$adminpc<max])
colnames(plotadmin)[1] <- "adminpc"

#pdf("effect_square.pdf",height=7, width=7*1.62)
ggplot(plots,aes(adminvec, meanvec))+stat_smooth(geom = "line")+
  theme_bw()+
  theme(axis.text=element_text(size=17),axis.title=element_text(size=17,face="bold"))+
  xlab("N. Top Tier Regional Govs, per 1m Citizens")+ylab("E(Expanded Services Index)")+geom_ribbon(aes(ymin=lowervec, ymax=highervec),alpha=0.2)+
  geom_point(data=plotadmin, aes(x=plotadmin$adminpc, y=min(lowervec)), alpha=0.2)+
  scale_x_continuous(breaks=c(0,10,20,30,40,50,60,70,80,90))
